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We study how the electronic structure of the bilayer graphene (BLG) is changed by electric field 
and strain from ab initio density-functional calculations using the LMTO and the LAPW methods. 
Both hexagonal and Bernal stacked structures are considered. The BLG is a zero-gap semiconductor 
like the isolated layer of graphene. We find that while strain alone does not produce a gap in the 
BLG, an electric field does so in the Bernal structure but not in the hexagonal structure. The 
topology of the bands leads to Dirac circles with linear dispersion in the case of the hexagonally 
stacked BLG due to the interpenetration of the Dirac cones, while for the Bernal stacking, the 
dispersion is quadratic. The size of the Dirac circle increases with the applied electric field, leading 
to an interesting way of controlling the Fermi surface. The external electric field is screened due to 
polarization charges between the layers, leading to a reduced size of the band gap and the Dirac 
circle. The screening is substantial in both cases and diverges for the Bernal structure for small 
fields as has been noted by earlier authors. As a biproduct of this work, we present the tight-binding 
parameters for the free-standing single layer graphene as obtained by fitting to the density-functional 
bands, both with and without the slope constraint for the Dirac cone. 

PACS numbers: 81.05.Uw; 73.22.-f 



I. INTRODUCTION 

Recently bilayer graphene (BLG) has been shown to 
possess a band gap in the presence of an external elec- 
tric field^*^ an effect that has important implications for 
transport across the graphene layers and for possible de- 
vice applications. While the freestanding BLG is a zero 
band gap semiconductor— like the single layer graphene 
(SLG), an appHed electric field through an external gate 
induces an asymmetric potential between the graphene 
layers, which in turn opens a gap between the valence 
and the conduction bands in the Bernal structure4ii^i^i^ 
It has been recently noted^" that graphene exhibits the 
hexagonal stacking surprisingly often, surprising because 
of the higher energy of the hexagonal structure. In the 
hexagonal structure, the two graphene layers are stacked 
vertically on top of one another, while in the Bernal str- 
cuture, one layer is rotated with respect to the other as 
indicated in Fig. [l] so that half of the atoms are directly 
above the carbon atoms in the first layer, while the re- 
maining half lie on top of the hexagon centers. The dif- 
ference in symmetry between the Bernal and the hexago- 
nal structures leads to substantial differences in the band 
structure, especially, in the formation of a band gap, an 
effect we study in this paper. 

To tune the field-induced band gap of the BLGs for 
practical applications in nano-electronic devices, it is im- 
portant to examine the factors having strong influence on 
its electronic structure. Electronic structure of BLG was 
shown to be sensitive to the interlayer spacing due to cou- 
pling between the graphene layer a^^i^^ . Hence, a uniaxial 
strain, which changes the interlayer spacing, could mod- 
ulate the field induced band gap ^^i^^ . The other factor 
that controls the gap is the screening of the electric field. 
An applied electric field causes unequal charge distribu- 
tion among the graphene layers which in turn creates a 




Bernal Stacking Hexagonal Stacking Interlayer Spacing (A) 

FIG. 1: Band structure of Bernal stacked (a) and hexago- 
nal stacked (b) BLG without any strain or electric field as 
obtained from the LMTO calculations. Fig. (c) shows the 
computed total energy per carbon atom as a function of the 
interlayer spacing. 



polarized electric field in the opposite directionji thereby 
reducing the effective electric field. A detail investigation 
through density-functional calculations is required to un- 
derstand the effect of the electric field and strain on the 
electronic structure. 

While the electronic structure of the Bernal BLG has 
been theoretically examined, such studies have not been 
performed for the hexagonal structure to our knowledge. 
In this paper, we report results from density-functional 
calculations and tight-binding models to examine the 
modification of electronic structure of both Bernal and 
hexagonal stacked BLG by electric field and strain. In 
the Bernal structure, the screening of the external electric 
field substantially reduces the band gap as shown earlier. 
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while in the hexagonal structure the screening is not as 
strong, as explained from a simple tight-binding model. 
In the Bernal structure, where a gap can be opened up by 
the electric field, strain can be used in addition to mod- 
ify its magnitude. We find that the gap can be increased 
considerably by reducing the interlayer spacing from its 
equilibrium value by a small amount. In the hexagonal 
structure, the gap cannot be opened up; however, an in- 
teresting point for this structure is that the topology of 
the bands leads to Dirac circles centered about the K 
point in the Brillouin zone with linear dispersion due to 
intcrpentration of Dirac cones. Though the electric field 
does not produce a gap, it increases the size of the Dirac 
circles leading to a novel way of controlling the Fermi 
surface. 

Density functional calculations presented in this pa- 
per were performed using the Linear MufHn-Tin Orbital 
(LMTO) Methodic and linear augmented plane wave 
(LAPW) method^^ using the local density approxima- 
tion (LDA). While the LAPW method was used to ob- 
tained the total energies and the optimized structures, 
the electronic structure was investigated using the LMTO 
method, which did not produce any substantial difference 
as far as the band structure is concerned. To study the 
field modulated band structure of the BLG, an extra po- 
tential A was added in the LMTO method to the on-site 
energies of the carbon s and p orbitals of one individual 
layer. We used the periodic boundary condition normal 
to the BLG planes, keeping about 10 A of vacuum be- 
tween the different BLG layers and, furthermore, in the 
LMTO method, we included layers of empty spheres at 
positions compatible with the symmetry of the particular 
BLG structure. The symmetry compatibility was espe- 
cially important for studying the screening effects at low 
electric fields. 



II. TIGHT-BINDING PARAMETERS FOR A 
SINGLE GRAPHENE MONOLAYER 

Before we present the results for the BLG, in this 
section, we present the tight-binding fitting parameters 
for the monolayer graphene obtained by fitting with the 
LAPW bands. These parameters will be used in a tight 
banding description of the BLG. 

For the freestanding monolayer graphene, it is known 
that the formation of the Dirac cones at the Fermi sur- 
face is the outcome of the Pz-Pz tt hopping between the 
two inequivalent carbon atoms A and B in the unit cell. 
Though a simple tight-binding model involving the Pz- 
Pz interaction within the nearest neighbor (NN) coordi- 
nation explains the linear band dispersion at the Dirac 
points K and K', it is essential to include further neighbor 
hoppings if an accurate description of the band structure 
in the entire Brillouin zone is desired. We find that for 
this, al least three NN hopping integrals must be retained 
as indicated from the root-mean-square deviation given 
in Table I and the band structure of Fig. ^ 



The TB Hamiltonian is a 2 x 2 matrix 
V h\s Hbb J 

where A and B denote the two carbon atoms in the unit 
cell. Retaining only the first two NN hoppings, the form 
of the matrix elements are 

hAA — <2[2cos(\/3fc2:a) -I- 4cos(\/3fc2;a/2) cos(3fcj,a/2)], 
hAB — ti[2exp(— i/cj^a/2) cos(\/3fc2;a/2) + exp(i/cj,a)], 

(2) 

and hAA = hsB- This can be easily generalized to further 
near neighbors. The symbols ti, t2, is, and ^4 are the NN 
hopping parameters as shown in Fig. [21 'a' is the C-C 
bond length, and the on-site energy of the carbon orbital 
is taken to be zero. The slope of the Dirac cone centered 
at the K or the K' point, easily obtained by diagonalizing 
the Hamiltonian and taking the limits, is given by the 
expression 

1; EE [dek/dk]K = a(3ii/2 - 3*3 + 9*4/4). (3) 

The slope at the Dirac point is about 5AeV- A, as cal- 
culated from LAPW, which corresponds to the velocity of 
8.2 X 10^ m/sec. The TB bands were fitted to the LAPW 
bands by keeping hopping integrals for a certain number 
of NNs and neglecting the hopping for further neighbors. 
We obtained two sets of optimized parameters. The first 
set of TB parameters was obatined by constraining the 
slope of the linear bands at the Dirac point to the LAPW 
value, while the second set was obtained without this con- 
straint. For low energy properties, where the magnitude 
of the slope may be important, the second set of the pa- 
rameters should be used. The results are shown in Fig. [2] 
and Table I. In our TB work of the BLG below, we have 
retained only the first NN interaction with the slope con- 
straint, since this is sufficient for our purpose as we are 
mainly interested in states close to the band gap region. 

III. BILAYER GRAPHENE IN THE 
HEXAGONAL STRUCTURE 

As mentioned already, recent experiments suggest that 
the hexagonal BLG is surprisingly common^ in spite of 
its higher energy as obtained from the total energy calcu- 
lation (Fig. [T]and Ref. 17) that indicate the Bernal BLG 
to be higher than the hexagonal BLG by ~ 5 meV/C 
atom. For the equilibrium interlayer spacing (d = 3.51 
A), the band structure for the hexaoganal stacked BLG 
is shown in Fig. El^a) and the topology of the bands near 
the Fermi surface is shown in Fig. [SI Unlike the case of 
the Bernal stacked BLG, here we see two interpenetrat- 
ing Dirac cones, which intersect to form a circular Fermi 
surface, the Dirac circle. 

An electric field does not change the overall feature 
of the band structure (see Fig. [3Kb)) and in contrast 
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FIG. 2: (color online) Tight-binding fitting (red and blue 
lines) of the LAPW bands (black lines) for a monolayer 
graphene and without the slope constraint at the Dirac point 
(see text). Tight-binding fitting was made for the p^ bands 
by retaining up to four nearest neighbors and the fitting pa- 
rameters are presented in Table 1. While just the first NN 
hopping is enough to fit the slope of the Dirac cone at the 
K point, at least three NNs must be retained to describe the 
band structure well in the entire Brillouin zone. 
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FIG. 3: Band structure of the hexagonal stacked BLG without 
(a) and with (b) an electric field and the variation of the 
radius of the Dirac circle (c). The charge density difference 
Sn between the layers induced by the electric field is shown 
in (c). All results were obtained from the density- functional 
LMTO calculations. 



TABLE I: Tight-binding NN hopping integrals obtained from 
the least square fitting of the LAPW bands. Fitting was done 
both with and without constraining the slope of the Dirac 
cone to its LAPW value. Quality of the fit is indicated by the 
root-mean-square deviation over the entire Brillouin zone. 
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FIG. 4: (color online) Topology of the band structures of 
hexagonal BLG with equilibrium interlayer separation and in 
the presence of an electric field. The Fermi surface is formed 
by two interpenetrating cones, forming a "Dirac" circle cen- 
tered around the K and the K' points in the Brillouin zone. 
The size of the circle can be controlled by the electric field as 
the right hand figure indicates for two diflferent electric fields. 



to the Bernal BLG, it docs not produce a band gap be- 
cause of the higher symmetry situation in the hexagonal 
stacking. The main features of the band structure can 
be understood from a tight-binding (TB) model involv- 
ing the nearest-neighbor tt- interaction of the Pz orbitals 
in the graphene plane and the a- interaction between the 
planes. With the four inequivalent carbon sites, one can 
form the Hamiltonian as in the monolayer case and in 
addition add the interlayer coupling term. One can then 
expand the Hamiltonian matrix elements for a small mo- 
mentum around the Dirac point, which then becomes an 
excellent approximation for the band gap region. The 



result is 

/ft i/^t \ 

where v = |tia is the Fermi velocity for the monolayer, 
t: — kx + iky is the complex momentum with respect 
to the Dirac point K, A is the potential difference be- 
tween the two layers, t is the interlayer hopping integral 
(~ 0.6eV as extracted from the LAPW results), and the 
basis set of the Hamiltonian consists of the Bloch func- 
tions made out of the carbon orbitals located at A, A, 



4 



B, and B, respectively and in that order. The atoms 
are indicated in Fig. [TJ The electric field potential A 
appearing in the Hamiltonian Eq. [4] corresponds to the 
final screened field experienced by the electrons and not 
to the bare electric field which is reduced by the polar- 
ization field due to screening, so that A — Aext/e, where 
Aext is the external electric field and e is the static di- 
electric constant. 

The Hamiltonian is easily diagonalized to yield the 
eigenvalues 



e{k) =±^±z/fc, 



where 



(5) 



(6) 



specifies the energy locations of the vertices of the Dirac 
cones. The energy dispersion is sketched in Fig. [6l which 
shows the two Dirac cones with linear dispersion with the 
vertices of the cones separated by the energy 2^. 

The radius of the circular Fermi surface, the "Dirac" 
circle, where the cones intersect the zero of energy, is sim- 
ply pq — £,/v. Since for small fields the screened potential 
A is proportional to the applied field, the equation indi- 
cates that po increases both with the applied electric field 
and the interlayer hopping parameter [t increases as the 
intcrlayer separation decreases). This is a notable fea- 
ture of the hexagonal BLG, viz., that we have a circular 
Fermi surface with zero band gap and that the radius 
of the circle can be modified by electric field and strain. 
Any doped holes or electrons will form a thin shell in the 
momentum space around the Dirac circle. This may lead 
to interesting electronic and magnetic properties in the 
presence of impurities. 



IV. BILAYER GRAPHENE IN THE BERNAL 
STRUCTURE 

The effect of strain and electric field on the band struc- 
ture of the Bernal BLG is summarized in Fig. [5] As 
seen from Fig. [Ha), for the Bernal BLG, there are 
four nearly parabolic bands near the Fermi energy and 
two of these touch at the Dirac point making it a zero 
band-gap semiconductor. However, in the presence of an 
electric field, these two bands split opening up a small 
gap (Fig. consistent with earlier density functional 

resultsi^iii^. The magnitude of the band gap increases 
linearly for small electric fields and saturates for higher 
electric fields as seen from Fig. [SJd). 

For large interlayer spacing, there is very little coupling 
between the two individual sheets and the bands for the 
single graphene sheet arc reproduced (Fig. ETd)). If the 
interlayer separation d is large enough (> 4.5 A), the gap 
vanishes irrespective of the electric field and decreasing 
the magnitude of d increases the gap substantially. Quite 
interestingly, if d is too small {d < 2.7 A), the band gap 
becomes indirect due to interactions with other orbitals 
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FIG. 5: Electric field modulated band structure of tlie Bernal 
stacked BLG for different interlayer spacings (a - c) and the 
variation of the band gap as a function of strain and electric 
field (d). Fig. (e) shows the magnitude of the polarization 
field Ep and the net screened electric field Es as a function 
of the applied electric field E. All results were obtained from 
the density-functional LMTO calculations. Topology of the 
bands in the gap region is sketched in Fig. (f). 



in addition to p^. However, the indirect band gap may 
not be realized in practice because of the large strain 
needed, which may also lead lead to a possible deforma- 
tion of the graphene sheets. For reasonable layer separa- 
tions, the electric field produces a small direct gap at k 
points, the locus of which is nearly a circle (the "Dirac" 
circle) centered around the K or the K' points in the Bril- 
louin zone as indicated from the band structure of Fig. 
Mh) and Fig. H 

The main features of the band structure can again be 
explained by a nearest-neighbor TB model, analogous 
to the case for the hexagonal structure. The resulting 
Hamiltonian is 
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(7) 



where the basis set of the Hamiltonian consists of the 
Bloch functions made out of the carbon orbitals located 
at A, B, A, and B, in that order. Diagonalization of the 
Hamiltonian yields the four eigenvalues 



e{k) = ± 



A^ + 4i/^7r7r'^ 



2r 



±2(4z.27r7rt(A2+i2)+t4)l/2]l/2^ 



(8) 



At the K or the K' points in the Brillouin zone, the en- 
ergies of the bands are e = ±^-^ -I- i^, and they 
disperse quadratically as seen from Fig. ^h) such that a 
gap of magnitude Eg — At/(A^ -I- 1^)^/^ is produced at 
k points in the Brillouin zone that make the "Dirac cir- 
cle" of radius po = A/2j/. The above expression clearly 
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shows that Eg increases linearly for low fields (A << t) 
and saturates at high fields (A >> t). This is consistent 
with the density-functional results for the dependence of 
the band gap on the external field E and bilayer spacing 
shown in Fig. [Hd). 



V. SCREENING IN THE BLG 

It has already been pointed out that the screening is 
divergent in the Bernal structure for small applied fields 
between the layers,**'^ which is a consequence of the band 
topology at low energies. In this Section, we examine the 
screening in the hexagonal structure and compare with 
the results for the Bernal structure. An applied electric 
field E transfers electrons from one graphene layer to the 
other, which produces a polarization field Ep resulting in 
the net screened field Eg that the electrons see, so that 
we have Eg = E — Ep. 

We compute the screening effects in three different 
ways using (a) density-functional LMTO method, (b) the 
tight-binding model, and finally (c) from an analytic ex- 
pression obtained by using the linear band dispersion. 
Unlike the Bernal structure, we find that the screening 
does not diverge in the hexagonal structure. 

In the LMTO calculations, an extra potential ±A/2 
was added to the on-site energies of the carbon orbitals 
of the two layers of BLG and the band structure was 
determined self-consistently. The screened potential A^ 
can then be directly obtained from the final on-site ener- 
gies of the carbon orbitals by examining the band-center 
energies. For this purpose, the C 2s orbital is especially 
useful, since it acts like a core orbital. 

For the TB as well as the analytic calculation, we com- 
puted the polarization field Ep by approximating the in- 
duced charges in the graphene layer by a uniform sheet 
charge density a, which turns out to be an excellent ap- 
proximation for computing the screening. The polariza- 
tion field is then given from the Gauss' Law in electro- 



statics to be Ep = cr/(2e), so that 



E, 



E — En 



E 



6n 
27' 



(9) 



where Sn is the difference between the sheet carrier den- 
sity of the two graphene layers and e is the dielectric 
constant. In terms of the electric potential the above 
equation can be wrtitten as 



A., = A - 



zdSr. 
"27 



(10) 



where d is the interlayer separation. One can obtain the 
value of Sn from the normalized eigenfunctions of the 
BLG Hamiltonian. 

Diagonalizing the Hamiltonian for the hexagonal BLG 
(Eq. [¥]), the eigenfunctions corresponding to the upper 
and the lower Dirac cones, centered at ^ and — ^ respec- 
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FIG. 6: The topology of the band structure of hexagonal 
BLG. The interpenetrating Dirac cones form a Dirac circle 
with radius po that makes the Fermi surface. The shaded 
region contributes to the integral in Eq. [T3]in the screening 
calculation. 



tively, are found to be 



p- 



( ±e- 

-9- 

V 1 



i0 \ 



I 



±9- 



V 1 

(1 ± A)i/2/2, 



where = tan '^{ky/k.^),p± = (1 ± ^)"'T^,q± 
(2t)^^(2^ ± A), and the ± sign inside the ket indicates 
the upper and the lower band and u{l) denotes the upper 
(lower) Dirac cone. Denoting the four components of the 
wave function by V'^, V'^, V'^, and ■0^ (from top to bot- 
tom), corresponding to the appropriate Bloch functions, 
the surface electron density is given by 



uk 



MA) 



B(B) 



(12) 



where Ace/; is the unit cell area, v is the band index, the 
superscript 1(2) refers to the two individual layers, and 
the factor of two in front of the summation takes care of 
the electron spin. 

The summation in Eq. [T^] need to be performed only 
over the shaded region in Fig. [51 since we are interested 
only in the difference — . The remaining portions of 
the two occupied bands ) and j-i/jL) cancel each other's 
contribution at each k point to this difference as can be 
easily verified from Eqs. pT|) and The net result is 
then 
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FIG. 7: The ratio between the apphed electric potential (A) 
and the screened electric potential (A3) as a function of A 
for the hexagonal and the Bernal stacked BLG. Results are 
shown from the analytic calculation, ( Eq. (|15|) for hexagonal 
and Eq. (|16|l for Bernal), the TB calculation, as well as from 
the density-functional LMTO calculation. 



where the summation goes over the two bands below Ep 
and the integration goes up to the radius of the Dirac cir- 
cle, so that the contribution comes only from the shaded 
region of Fig. [HI The factor of four in Eq. [T3] was in- 
serted to account for both the spin degeneracy and the 
two Dirac valleys in the Brillouin zone at the K and K' 
points. Using the wave functions Eq. [Til the integration 
is easily performed to yield 



dn 



2^ 



(14) 



Now, electrons of course see the screened field, so that 
the net charge density difference 5n between the layers is 
obtained by using the screened potential in lieu of A 
in the expression ITU) . This together with Eq. [TU] yields 
the ratio between the external potential and the screened 
potential, which is 



A_ 

a: 



1 



Airev 



V4i2 + A2. 



(15) 



This is plotted in Fig. [H] together with the LMTO 
result and the TB result. The latter was obtained by a 
full self-consistent calculation using the nearest-neighbor 
tight-binding model, from which the charge difference dn 
was computed and used in Eq. [10] to obtain the ratio of 
A/ As- In the LMTO results, all electrons including the 
carbon a electrons participated in the screening, while 
in the analytical and TB calculation, which included the 
effects of only the tt electrons, we replaced, for simplicity, 
just the vacuum dielectric constant for e appearing in Eq. 

m 



Fig. [Tjshows the calculated results for the screening for 
both the hexagonal and the Bernal structures. For both 
cases, screening increases with applied electric field, ex- 
cept for the divergence at low field for the Bernal struc- 
ture as has been pointed out earlier. The increase of 
screening with increasing field is due to the fact that 
proportionately more carriers become involved in screen- 
ing as the density-of-states for the linear bands near the 
Fermi energy increase as square of energy. The analytic 
and the TB results agree quite well, which is not surpris- 
ing since the energy spectrum differs only beyond the lin- 
ear region of the bands and these states don't contribute 
much because they are far away from the Fermi energy. 
The density-functional LMTO result, on the other hand, 
shows larger screening because it includes all carbon elec- 
trons that participate in screening. 

For the Bernal structure, the divergent screening for 
low fields may be explained due to the relatively flat 
bands at the Fermi energy (see Fig. [T]a). The effects of 
the two flat bands touching at the Fermi energy may be 
taken into account by constructing a 2 x 2 effective Hamil- 
tonian using perturbation theory from the full Hamilto- 
nian of Eq. [7| and then computing the magnitude of the 
charge difference from the eigenfuctions. This can be 
done analytically and, for low fields, the result is^ 



A 

a: 



2dm , 2fc^ 
— TT- In 



mA, ' 



(16) 



where m = ti(2v'^)~^ and the cutoff momentum kc ~ 
0.065 A~^. This result has been plotted in Fig. [Tjfor the 
Bernal structure and it explains the calculated trend in 
the density-functional results. 



VI. SUMMARY 

In summary, we have studied the effect of the electric 
field and strain on the electronic structure of hexagonal 
and Bernal stacked bilayer graphene. While a gap opens 
up in the Bernal structure BLG by the application of 
an external electric field, the hexagonal BLG remains a 
zero-gap metal, but with an interesting circular Fermi 
surface, the Dirac circle, that can be controlled by both 
strain and electric field. Any doped electrons or holes 
will occupy a thin shell in the Brillouin zone, which leads 
to an interesting Fermi surface in the doped system. 

This work was supported by the U. S. Department of 
Energy through Grant No. DE-FG02-00ER45818. 
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